Effects of edge disorder on the stability of quantum oscillations in two-dimensional coupled systems

This paper utilizes the theory of quantum diffusion to analyze the electron probability and spreading width of a wavepacket on each layer in a two-dimensional (2D) coupled system with edge disorder, aiming to clarify the effects of edge disorder on the stability of the electron periodic oscillations in 2D coupled systems. Using coupled 2D square lattices with edge disorder as an example, we show that, the electron probability and wavepacket spreading width exhibit periodic oscillations and damped oscillations, respectively, before and after the wavepacket reaches the boundary. Furthermore, these electron oscillations exhibit strong resistance against disorder perturbation with a longer decay time in the regime of large disorder, due to the combined influences of ordered and disordered site energies in the central and edge regions. Finally, we numerically verified the universality of the results through bilayer graphene, demonstrating that this anomalous quantum oscillatory behavior is independent of lattice geometry. Our findings are helpful in designing relevant quantum devices and understanding the influence of edge disorder on the stability of electron periodic oscillations in 2D coupled systems.

Quantum dynamics 30,31 are crucial for understanding the electronic transport properties in ordered and disordered systems.In quasi-1D and quasi-2D systems with Anderson disorder, although the electronic wave functions cannot extend to the boundaries, they can still hop regularly between chains or layers.In our previous research 29,30 , we used the method of quantum dynamics to predict that in low-dimensional coupled systems with mirror symmetry, periodic electronic oscillations exist.When all site energies are randomly disordered, these periodic electronic oscillations turn into damped oscillations with a continuously decreasing decay time as the degree of disorder increases.However, none of these studies have addressed the effects of edge-site energy disorder on the stability of quantum oscillations.Ultrathin quantum films often exhibit roughness, vacancies, and defects at edges due to cutting, material growth, adsorption, and other factors during materials fabrication.The random disorder of energies at the edge sites could affect the electronic transport properties.Therefore, it is interesting to study the stability of electron oscillations in 2D coupled systems with edge-site energy disorder.
This research is based on the Anderson model of 2D coupled systems.It employs the tight-binding approximation theory and the method of direct diagonalization of Hamiltonian matrix to solve both the static Schrödinger equation and the time-dependent Schrödinger equation.Additionally, the quantum diffusion theory is utilized to derive the time evolution of electron probability and wavepacket propagation width on each layer.We use a coupled bilayer of square lattices to explore the effects of edge-site energy disorder on the stability of the periodic electron oscillations in 2D coupled systems, and use the commonly studied material, bilayer graphene, to confirm the universality of these conclusions.

Theory
The proposed model is a finite-sized 2D coupled square lattices, where the edge site energies follow Anderson disorder, as illustrated in Fig. 1.During the numerical simulations, we use the single-electron tight-binding approximation theory, employing the method of direct matrix diagonalization and fixed boundary conditions to analyze the electron diffusion behavior in finite-sized 2D coupled systems.The Hamiltonian in the tight-binding form for a system with 2N sites is given by where H 1 , H 2 represent the sub-Hamiltonians of the layer 1, layer 2, and H int represents the interlayer Hamilto- nian between the two-layers.Based on the tight-binding approximation theory, H 1 , H 2 and H int can be described as follows where |1, i� , |2, i� represent the basis sets, and ε i represent the onsite energies at the i-th site in layer 1 and layer 2, respectively.This paper only considers the effects of disorder in edge-site energies, so the energies of edge-sites (colored in Fig. 1, denoted as ε out ) are taken from random numbers in the interval [− W,W ], referred to as the disordered region, while the energies of sites in the center region (gray in Fig. 1, denoted as ε int ) are set to 0, referred to as the ordered region, where W represents the degree of disorder and corresponds to the periodic system with W = 0. To simplify the calculation, only the hopping parameters of nearest-neighbor sites for intralayer hopping , where i and j represent the i -and j-th sites.During the numerical simulations, all physical (1) quantities are expressed in dimensionless units.The lattice constant and the intralayer hopping strength are normalized with a = 1, h = 1 and u = 1 , respectively.The interlayer hopping strength u and the degree of disorder W are treated as modulated variables for tuning quantum diffusion.
The electron eigen-states of the 2D coupled system can be obtained by the following static Schrödinger equation where E and | ψ� are the eigen-energy and eigen-wavefunction, respectively.If layer 1 and layer 2 have a mirror- symmetric structure, regardless of edge-disordered or not, H 1 and H 2 are iso-spectral with the same spectrum of electron eigen-energy e l (l = 1, 2, . . .N) , and their corresponding electron wavefunctions |ϕ (1) l , |ϕ (2) l are respectively given by Projecting | ψ� to the basis set formed by |ϕ (1) l and |ϕ (2) l , one has Substituting Eq. ( 7) into Eq.( 5) and using Eqs.(1)-( 4), one can derive and describe the electron wavefunctions and eigen-energies of the 2D coupled system as follows The dynamics of an electron wavepacket in the 2D coupled system are governed by the time-dependent Schrödinger equation where |� (t)� = (� 1 , � 2 , . . .� 2N ) T represents the electron wavefunction at time t , with � n (t) being the coef- ficient of |�(t) � corresponding to site n .If the initial state of electron wavepacket is at position m in layer 1, then |� (t)� t=0 = |1, m� , and thus where, C * lm is the complex conjugate of C lm .Projection of |�(t) � to the basis set of layer 1 and layer 2 leads to where � (s) is the wavefunction of the single isolated layer with initial distribution of wavepacket 31 .The probability of finding electron at the n-th site on layer 1 and layer 2 follows where is the probability of finding electron at site n on the single isolated layer.As the probability of finding electron on layer 1 is given by P 1 (t) = N n=1 P (1) n (t) , on layer 2 is given by P 2 (t) = N n=1 P (2) n (t) , and taking into account the rule of P 1 (t) + P 2 (t) = 1 , one has From Eqs. (13a) and (13b), one can also derive that the spreading widths of electron wavepacket d 1 (t) on layer 1, defined as d 2 1 (t) = N n=1 P (1)  www.nature.com/scientificreports/with where d 2 nm (t) is the spreading width of wavepacket on isolated single layer without coupling, and n , m represent the planar positions of the n-th and m-th sites, respectively.The spreading width of electron wavepacket of the 2D coupled system is described by According to the theory of quantum diffusion, the relationship between d(t) and time t follows d(t) ∝ t β with 0 ≤ β ≤ 1 , where β = 0, 0.5, 1.0 corresponds to localization, normal diffusion, and ballistic diffusion, respectively, 0 < β < 0.5 and 0.5 < β < 1, correspond to sub- diffusion and super-diffusion.The analysis shows that, only in the special case of strong localization does the localization length become close to the lattice spacing, d nm (t) ∼ a .At this moment, d 1 (t) and d 2 (t) exhibit dif- ferent amplitudes of oscillation.However, in general, d nm (t) is much larger than a , which means that, for large time d nm (t) ≈ d(t).In this case, the spreading width of wavepacket or the oscillation length on each layer can be given by By defining the spreading width of wavepacket ratio d k (t) = d k (t)/d(t) one can obtain Equations (14a), (14b), and Eqs.(18a), (18b) clearly indicate that both the probability P k (t) of finding elec- tron on k-th layer, and the square of the spreading width of wavepacket ratio d 2 k (t) , exhibit a standard periodic oscillation with frequency f = u/π , which is determined only by the interlayer hopping strength.This means that the 2D coupled system with a mirror-symmetric structure, regardless of order or edge-disorder, exhibits a periodic oscillation in both transverse and longitudinal directions with the same frequency, due to the iso-spectra of two-layers, which are independent of system size and type of intralayer potentials.

Numerical results
The numerical results for the edge-disordered bilayer square lattices and the bilayer graphene are as follows.In the numerical calculations, the initial position of the electron wavepacket is located at the central site of layer 1. Considering both the accuracy of conclusions and the statistical significance, the results presented below are from the average of 20 different random samples, each with a system size of 2 N = 10,082 for the bilayer square lattices and 2 N = 7688 for the bilayer graphene.Since P 1 (t) + P 2 (t) = 1 and d 2 1 (t) + d 2 2 (t) = 1 , we solely focus on the analysis of P 1 (t) and d 2 1 (t).Figure 2 illustrates the electron diffusion in an edge-disordered coupled system of bilayer square lattice (BSL).In Fig. 2a, it can be observed that the spreading width of wavepacket, denoted as d(t) , exhibits linear increase with time before the wavepacket reaches the boundary (t < 15.7), and the increase ceases after the wavepacket reaches the boundary (t > 15.7), which is similar to the case without edge disorder 30 .The dashed line in the P b (t) inset marks the time when the wavepacket reaches the boundary, and the time intervals 0 ≤ t ≤ 10 and 75 ≤ t ≤ 85 of the d k (t) and P k (t ) insets correspond to the scenarios before and after the wavepacket reaches the boundary.From the figures of d k (t) and P k (t ) sets one can see that the amplitude of the electron oscillations significantly decreases after the wavepacket reaches the boundary.In Fig. 2b, we present the Fast Fourier transformation (FFT) spectrum of the electron probability on layer 1, P 1 (t) , under different disorder strength W .The figure shows that, as W increases, P 1 (t) maintains a single-frequency signal pattern, and the signal strength undergoes a transition from a decrease to an increase at W = 3, marking a difference from traditional Anderson disordered systems 30 .
Numerical analysis reveals that in the edge-disordered 2D coupled system, before the wavepacket reaches the boundary, both the electron probability P k (t) and the square of the spreading width of the wavepacket ratio d 2 k (t) exhibit periodic oscillations, similar to those observed in periodic systems, following the Eqs. ( 14 and (18).After the wavepacket reaches the boundary, both P k (t) and d 2 k (t) display damped periodic oscillations similar to traditional Anderson disordered systems, as described by the following where, t 0 represents the decay time, V /π is the frequency parameter, and t b is the time when the peak of the wavepacket reaches the boundary.www.nature.com/scientificreports/ Figure 3 illustrates the electron oscillations of probability P 1 (t) and the square of wavepacket spreading width ratio d 2 1 (t) in an edge-disordered bilayer square lattice.Figure 3a and b indicate that, before the wavepacket reaches the boundary ( t ≤ 15.7), both P 1 (t) and d 2 1 (t) exhibit periodic oscillations, following Eq. (14a) and Eq.(18a), respectively.After the wavepacket reaches the boundary ( t > 15.7), these peaks of periodic oscillations undergo exponential decay over time and become damped oscillations, which can be described as Eq. ( 19).More importantly, the decay rate 1/t 0 of these damped oscillations exhibits a transition from acceleration to deceleration as W increases, occurring at W = 3.In the two figures, W = 1.2, 1.8, 3.0, 8.0 correspond to t 0 = 99.4,44.2, 21.5, 48.4 and V = 1.00171, 1.00289, 1.00047, 0.99779 , respectively.Here, t b = 15.7 which can be inferred from P b (t) .The fitting parameters t 0 , t b , V for d 2 1 (t) are the same as those for P 1 (t) .According to Eq. ( 19), one can infer that ln(2P 1 − 1) = t/t 0 .Figure 3c presents the numerical results of ln(2P 1 − 1) as a function of t for different W , and their corresponding linear fits of the linear equation ln(2P 1 − 1) = t/t 0 , where W = 0.6, 0.8, 1.0, 1.2, 3.0, 5.0, 8.0, 10 correspond to t 0 = 397.84,226.66, 144.19, 99.45, 21.46, 29.18, 48.43, 62.76.To obtain precise fitting results, we initially estimate the value of t 0 by performing a test fit using numerical data with large time values ( t > t 0 ).Subsequently, we extract the refined fitting parameters t 0 by conducting repeated fits using data dur- ing 0 < t < t 0 , until the slope 1/t 0 remains constant.In this figure we can see a significant reversal in the fitting slope at W = 3.As W increases, the slope increases for W < 3 and decreases for W > 3.By substituting t 0 into Eq.( 19), the frequency parameter V can be determined by fitting the data for 0 < t < t 0 of P 1 (t) .Figure 3d presents the numerical results of t 0 with different W in log-log and their linear fitting with function logt 0 ∼ log W , which clearly indicate that t 0 ∼ W α , with α = -1.85 for W < 3, and α = 1.07 for W ≥ 3, where t 0 is the decay time of the electron oscillation.This anomalous quantum diffusion phenomenon suggests that, in edge-disordered 2D coupled systems, the damped oscillations will improve with the degree of disorder increases when W exceeds a critical value, larger disorder leads to a longer decay time.The inset shows that the frequency V /π only undergoes a slight change with different degrees of disorder, V = V − V 0 ≈ 0 .All of those indicate that, compared to traditional disordered 2D coupled systems 30 , the electron oscillations in edge-disordered 2D coupled systems are more stable against disorder.
The anomalous quantum diffusion behavior in edge-disordered coupled systems results from the combined influence of disordered edge-site energies and ordered site energies in the central regions.The degree of electron localization in quantum diffusion can be described by the particle participation number, which is defined as , where φ(E, r n ) is the wave function of the state E .Here, P(E) is related to the system size 2N as P(E) ∝ (2N) γ , with γ = 1 for extended states, γ = 0 for localized states, and 0 < γ < 1 for critical states.Figure 4 shows that, the energy spectrum of the edge-disordered coupled system is divided into a stable central region and two extended localized tails, separated by the critical energy levels E c = ± 5, due to the influence of the ordered central part and the disordered edges.Figure 4a displays the density of states (DOS) for different values of W . Figure 4a shows that the DOS exhibits the same behavior as a bilayer periodic system for any values of W in the band center region (|E|< 5), while its band tails ( |E|> 5) continuously expand as W increases, similar to those in a traditional disordered bilayer 27 .We take the edge-disordered system with W = 5 as an example to demonstrate the relationship between P(E) and the system size 2N , as shown in Fig. 4b and c.From Fig. 4b, it is evident that the numerical values of P(E) noticeably increase as 2N increases for |E|< 5, but independent of 2N for |E|> 5.The detailed analysis in Fig. 4c shows that, there is a significant mobility edge at |E|= ± 5, dividing the energy spectrum into a ballistic central region ( |E|< 5, γ = 1) and two localized tails ( |E|> 5, γ = 0).Considering the symmetry of the P(E) spectrum and for the sake of graphical clarity, the figure only shows several typical energy levels for E > 0.
Consider a bilayer square lattice system with 2N sites.By solving the Schrödinger equation, the electron wavefunction component at the n-th site at time t can be given by www.nature.com/scientificreports/ Here, E j j = 1, 2, 3, 4, . . .2N are the eigen-energies, φ n E j is the eigen-wavefunctions corresponding to E j at n-th site, where φ * 0 E j is the complex conjugate of φ n E j .Therefore, the probability of finding electron at n-th  19) with red lines(t > 15.7), and both Eq.(14a) and Eq.(18a) with blue line(t < 15.7).(c) The numerical results of ln (2P 1 − 1)(symbols) as a function of t , and the corresponding fitting results (dash lines) with linear equation ln(2P 1 − 1) = t/t 0 , where P 1 corresponds to the amplitude of oscillations of P 1 (t) .(d) Values of log t 0 versus log W (symbols) and the corresponding linear fitting (dashed line) with log t 0 ~ log W , as well as values of V versus W (inset), where V = V − V 0 , and V 0 = u is the parameter of the bilayer square lattice without disorder.www.nature.com/scientificreports/site at time t can be expressed as According to the previous discussion, the electron wavefunc- tion in an edge-disordered bilayer system encompasses both extended states and localized states.Therefore, the probability of finding electron at k-th (k = 1, 2) layer can be expressed as Here, NO P no (t) represents the probability of finding electron in the central region with ordered sites ( ε int = 0) in the k-th layer, and ND P (k) nd (t) represents the probability of finding electron at the edge disordered sites ( ε out ∈ [−W, W]) in the k-th layer.According to the theory of electron localization, as W increases, the edge wavefunctions become more and more localized.In the case of large disorder, ND P (k) nd (t) approaches to a constant, and thus P k (t) → NO P (k) no (t) + constant, which means that the electron oscillation behavior will gradually approaches that in a bilayer periodic system.This is also the reason why, in the edge-disordered coupled system, when W exceeds a critical value, larger disorder leads to longer decay times in damped oscillations.This anomalous behavior is similar to the crossover of quantum diffusion in the bilayer system with order-disorder separation [26][27][28]32 .
However, compared to the square lattices, hexagonal lattices are more prevalent in 2D materials.In the following sections, we will take a bilayer graphene of AA-stacking, shown in Fig. 5, as an example to validate the universality of electron oscillations in edge-disordered 2D coupled systems.The system size of the bilayer graphene used for numerical computation is 2N = 7688 , with x = 31 × 2a , y = 31 × √ 3 2 a , z = d , where a is the lattice spacing equal to 0.142 nm 33 , and d is the distance between the layers equal to 0.335 nm 33 .Additionally, the nearest-neighbor hopping energies h is −2.7 eV 33,34 .To simplify the numerical calculation, all physical param- eters mentioned above are used in dimensionless form.Here, a , d , and h are set to 0.142, 0.335, -2.7, respectively.The nearest-neighbor interlayer hopping parameter u is set to 1, the site energies of the central ordered region are set to ε int = 0, while energies of the disordered sites at the edge are set to ε out ∈[−W, W ]. The numerical results presented below are from average of 20 different disordered configurations to ensure reliable results.The initial location of the electron is at the center site of layer 1.
The time evolution of the wavepacket in an edge-disordered bilayer graphene, as shown in Fig. 6, exhibits similar electron oscillations to those in bilayer square lattice.These oscillations can be described by Eqs. ( 14) and ( 18) before the wavepacket reaches the boundary, and by Eq. ( 19) after it reaches the boundary.Figure 6a illustrates the spreading width of the wavepacket d(t) with W = 3, which first exhibits linear growth with time and then ceases to grow after the wavepacket reaches the boundary.The insets show the time evolution of the probability of finding electrons at the edge sites P b (t) , indicating the time when the wavepacket reaches the boundary at t ≈ 5 , as well as the spreading width d k (t) and electron probability P k (t) on the k-th layer, which indicate the decrease in the amplitude of the electron oscillations over time.The FFT frequency spectrum of P 1 (t) , as shown in Fig. 6b, exhibits single-frequency signals for any degree of disorder W , and the signal strength undergoes a transition from decrease to increase at W = 8. Figure 6c and d exhibit the oscillations of P 1 (t) and d 2 1 (t) of typical disorder and their excellent fits.The fitting parameters for W = 1.8, 3.0, 8.0, 19 are t 0 = 71.115,36.112,21.858 and V = 1.00288, 1.00531, 1.00823 .The fitting parameter t 0 is obtained from the linear fitting of the damped oscillation amplitudes ( t > 6.28), as shown in Fig. 6e, with the fit function ln(2P 1 − 1) = t/t 0 .By substituting t 0 into Eq.( 19) and then fitting the damped oscillation data of P 1 (t) and d 2 1 (t) , we can get the parameter V. From Fig. 6e we can see a significant reversal in the fitting slope at W = 8.The relationship of t 0 and V versus W are shown in Fig. 6f, which clearly indicate that t 0 ∼ W α , with α = − 1.78 for W < 8, and α = 0.62 for W ≥ 8.The inset shows the parameter V only undergoes a slight change with different degrees of disorder.

Conclusions and discussions
In summary, this study finds that, in the coupled bilayer square lattices of edge disorder, the probability and spreading width of an electron wavepacket initially located on one layer exhibit periodic oscillations before the wavepacket reaches the boundary and become damped oscillations after the wavepacket reaches the boundary.As the disorder strength increases, the frequencies of these electronic oscillations remain almost unchanged, and  the decay time of damped oscillations exhibits a transition from a decreasing behavior to an increasing behavior of decay time at a critical disorder value, which implies that the electron oscillations exhibit strong resistance against disorder perturbation in the regime of large disorder.We confirmed the universality of this anomalous damped oscillation behavior using bilayer graphene systems, emphasizing its independence from the geometric shape of the lattice.We note that, although damped periodic oscillations have been reported in the case of internal disorder in previous work 29,30 , it is unknown whether the oscillations under edge-disordered bilayer systems have the same form of damped oscillations.The unique aspect of this paper is that our numerical results reveal that the oscillations maintain the same form of damped oscillations under both small and large edge disorder.Particularly, there exists a crossover from a worse damped oscillation to an improved damped oscillation as W increases, which is a new and unconventional phenomenon unreported in the literature.Since the site energies of edges are of the Anderson type of disorder given by [−W, W] , the edges have both infinity walls and infinity wells, when W goes infinity.Classically, some particles will be trapped in the wells, and some will bounce back and eventually will be trapped in the wells.Therefore, the phenomenon of the improved damped oscillations under large disorder is unexpected and anomalous from the classical point of view.Our results are significant in understanding the effects of edge-disorder on the stability of periodic electronic oscillations, which is important in designing novel quantum devices.

Figure 1 .
Figure 1.Schematic illustration of the edge-disordered 2D coupled system.h and u correspond to electron hopping parameters of intralayer and interlayer between nearest-neighbor sites, respectively.ε int and ε out represent the onsite energies in the center region(gray) and edge region(colored), respectively.The different colors indicate different types of atoms in the figure.

Figure 2 .
Figure 2. Electron diffusion in an edge-disordered bilayer square lattice with 2N = 71 × 71 × 2 , h = 1, u = 1,ε int = 0, ε out ∈ [−W, W ]. (a)The spreading width of electron wavepacket d(t) in the edge-disordered 2D coupled system with W = 3.The insets show the probability of finding electron at the edge-sites P b (t) , the spreading width d k (t) and the probability P k (t) of find electron in the k-th layer, where k = 1,2 corresponds to the blue and red curve, respectively.(b) FFT frequency spectrum of probability P 1 (t) as a function of W.

Figure 3 . 2 1 2 1
Figure 3. Electron oscillations of probability P 1 (t) and the square of wavepacket spreading width ratio d 2 1 (t) in an edge-disordered bilayer square lattice, with 2N = 71 × 71 × 2 , h = 1, u = 1, ε int = 0, ε out ∈ [−W, W ]. (a) and (b) Show the numerical results (open circles) of P 1 (t) and d 2 1 (t) for W = 1.2, 1.8, 3.0, 8.0 , and the corresponding fitting results based on Eq. (19) with red lines(t > 15.7), and both Eq.(14a) and Eq.(18a) with blue line(t < 15.7).(c) The numerical results of ln (2P 1 − 1)(symbols) as a function of t , and the corresponding fitting results (dash lines) with linear equation ln(2P 1 − 1) = t/t 0 , where P 1 corresponds to the amplitude of oscillations of P 1 (t) .(d) Values of log t 0 versus log W (symbols) and the corresponding linear fitting (dashed line) with log t 0 ~ log W , as well as values of V versus W (inset), where V = V − V 0 , and V 0 = u is the parameter of the bilayer square lattice without disorder.

Figure 5 .
Figure 5. Schematic illustration of bilayer graphene of AA-stacking(the upper atom is directly above the lower atom).Parameters h and u correspond to nearest-neighbor hopping energy of intralayer and interlayer, respectively.ε int and ε out are the site energies in the central ordered region and in the edge disordered region, respectively.

Figure 6 . 2 1
Figure 6.The electron oscillations in the edge-disordered bilayer graphene with 2N = 7688, h = −2.7,u = 1 , ε int = 0, ε out ∈ [−W, W] .(a) The spreading width of the wavepacket d(t) at W = 3.The insets show the probability of finding electron at edge-sites P b (t) , as well as the spreading width d k (t) and electron probability P k (t) in k-th layer, where k = 1 and 2 correspond to the blue and red lines, respectively. (b) FFT frequency spectrum of probability P 1 (t) as a function of W . (c) and (d) represent P 1 (t) and d 2 1 (t) for W = 1.8, 3.0, 8.0, 19 (open circles) and the corresponding fitting with Eqs.(14a), (18a) (blue lines) and Eq.(19)(red lines).(e) The numerical results of ln (2P 1 − 1)(symbols) and their corresponding linear fitting (dash lines) for different degrees of disorder W. (d) The values of t 0 versus W in log-log(symbols) and their corresponding linear fitting (dashed lines), indicating t 0 ∼ W α , with α = − 1.78 for W < 8, and α = 0.62 for W ≥ 8.The inset shows the values of V versus W, where V = V − V 0 .